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We address the problem of directional mobility of discrete solitons in two-dimensional rectangular 
lattices, in the framework of a discrete nonlinear Schrodinger model with saturable on-site non- 
linearity. A numerical constrained Newton-Raphson method is used to calculate two-dimensional 
Peierls-Nabarro energy surfaces, which describe a pseudopotential landscape for the slow mobility 
of coherent localized excitations, corresponding to continuous phase-space trajectories passing close 
to stationary modes. Investigating the two-parameter space of the model through independent vari- 
ations of the nonlinearity constant and the power, we show how parameter regimes and directions 
of good mobility are connected to existence of smooth surfaces connecting the stationary states. In 
particular, directions where solutions can move with minimum radiation can be predicted from flat- 
ter parts of the surfaces. For such mobile solutions, slight perturbations in the transverse direction 
yield additional transverse oscillations with frequencies determined by the curvature of the energy 
surfaces, and with amplitudes that for certain velocities may grow rapidly. We also describe how 
the mobility properties and surface topologies are affected by inclusion of weak lattice anisotropy. 

PACS numbers: 42.65.Wi, 63.20.Pw, 63.20.Ry, 05.45.Yv 



I. INTRODUCTION 

Intrinsically localized modes, or discrete solitons 
(breathers) , appear as generic excitations in a large vari- 
ety of physical systems [1 j, where spatial periodicity (or 
discreteness) provides gaps in the linear dispersion rela- 
tion and nonlinearity allows for detuning the oscillation 
frequencies into these gaps, resulting in spatial localiza- 
tion. Particularly important applications of discrete non- 
linear systems are nonlinear optical waveguide arrays [2] , 
with possibilities to change and control all essential pa- 
rameters, such as geometry, dimensionality, nonlinearity, 
beam angle, etc. For weakly coupled waveguides, the 
relevant mathematical description, derived via coupled 
mode theory |2j [3] , is a Hamiltonian lattice equation of 
the discrete nonlinear Schrodinger (DNLS) type [TH3], 
where the particular type of nonlinearity depends on the 
nonlinear response of the medium. The most studied case 
is a DNLS equation with cubic on-site nonlinearity [4], 
corresponding to a Kerr medium [2], which also appears 
generically as a modulational equation for the small- 
amplitude dynamics in chains of coupled anharmonic os- 
cillators pp. However, photorefractive media also enables 
the generation of discrete spatial solitons [2 , and the 
corresponding lattice model is a DNLS equation with 
saturable nonlinearity (s-DNLS) which was studied in a 
number of theoretical works [5HTT]. 

A particularly interesting property of the one- 
dimensional (ID) s-DNLS model [5]-[9] is the existence of 
certain "sliding velocities" , where localized discrete soli- 
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tons may travel in the lattice without radiation. This 
behavior was connected to the existence of "transpar- 
ent points" associated with the vanishing of a so-called 
Peierls-Nabarro (PN) potential barrier [12], usually de- 
fined as the difference in energy (Hamiltonian) at con- 
stant power (norm) between the two fundamental local- 
ized stationary solutions centered at one site (odd mode) 
and symmetrically in between two sites (even mode). 
Close to the points of vanishing of this energy difference 
are regions of stability exchange between the even and 
odd solutions, associated in general with the appearance 
of a family of intermediate, asymmetric stationary solu- 
tions, connecting both types of symmetric solutions at 
the bifurcation points [13 . The existence of regimes of 
enhanced mobility close to such bifurcation points is also 
well-known from studies of other one-dimensional lattice 
models [T3HT5] . 

On the other hand, it is still an open issue whether 
moving discrete solitons may exist as localized, radiation- 
less modes also in two-dimensional (2D) lattices. As was 
shown numerically in [10 , a scenario with exchange of 
stability through bifurcations with asymmetric station- 
ary solutions appears also for the 2D saturable model in a 
square (isotropic) lattice, involving in this case three dif- 
ferent types of fundamental solutions [16] : one-site (odd- 
odd, 00), two-site (odd-even, OE), and four-site (even- 
even, EE) modes. It was also shown numerically in [10] . 
that solutions with good (but generally not radiation- 
less) mobility in the axial directions may exist in these 
regimes, and that the necessary energy needed for render- 
ing a mobile stable stationary solution agreed well with 
the concept of PN-barrier, if its definition was extended 
to take into account also the energy for the relevant in- 
termediate stationary solution (an analogous situation is 
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well-known for the PN potential of kinks [17]). Very sim- 
ilar observations were also made later for a 2D DNLS 
model with cubic-quintic nonlinearity [18] . Moreover, in 
the low-power regime of the 2D s-DNLS equations, good 
mobility was also observed in diagonal directions [TP] , 

However, the relation between existence of regions of 
stability exchange, small PN-barrier and mobile localized 
solutions in 2D is not that trivial, as shown, e.g., for a 
model with cubic inter-site nonlinearities in [19] : even in 
regimes with small PN barrier and existence of interme- 
diate solutions, the mobility may be very poor, if there 
is no continuous path in phase space passing close to the 
relevant stationary modes. On the other hand, as was ob- 
served for the low-power (i.e., close to continuum limit) 
regime of a 2D lattice with quadratic nonlinearity in [20] , 
the effective Peierls-Nabarro potential may in some sit- 
uations be weak enough to allow mobility in arbitrary 
directions, without any direct connection to bifurcations 
and symmetry-broken stationary solutions. 

So there is clearly need for a better understanding 
of the conditions for mobility in 2D lattices. It is the 
purpose of the present paper to generalize the concept 
of PN-barrier as discussed above, and introduce a full 
2D PN potential surface describing the pseudo-potential 
landscape in-between all stationary modes. We will use 
a numerical constrained Newton-Raphson (NR) method, 
previously applied to ID lattices in [21J, to explicitly con- 
struct these surfaces for the 2D saturable model from [10] , 
and show how parameter regimes and directions of good 
mobility may be immediately identified from smooth, flat 
parts of these surfaces. We will also illustrate how the 
interplay between translational motion in one lattice di- 
rection and oscillatory motion in the orthogonal direction 
can be intuitively understood from the topology of the 
corresponding PN surfaces. 

The structure of the paper is as follows. In Sec. [TTJ 
we describe the 2D s-DNLS model, its stationary local- 
ized solutions (Sec. II A), and the explicit implementa- 



tion of the numerical constraint method used to calcu- 



late the PN potential surfaces (Sec. II B). We point out 
some essential differences between the 2D implementa- 
tion and earlier implementations for the ID case. In Sec. 
Ill A we identify the regimes of parameter space where 
smooth surfaces can be found for the isotropic lattice, 
and make a classification of appearing surfaces with dif- 
ferent topologies. In Sec. IIIB we illustrate with direct 



dynamical simulations various types of behavior for dis- 
crete solitons moving slowly in different directions, corre- 
sponding to the surfaces discussed in Sec. |III A| We also 
here show examples on how simultaneous excitation of 
transverse oscillation modes may affect discrete solitons 
moving slowly in axial directions. In Sec. |IV[ we ana- 
lyze effects of adding a weak anisotropy to the lattices 
in the parameter regimes used in previous sections, in 
order mainly to investigate whether this may promote 
mobility in directions different from the axial or diag- 
onal ones (which presumably would be associated with 
the appearance of additional "valleys" in the PN poten- 



tial surface). Finally, in Sec. [V] we make some concluding 
remarks, summarizing our results and pointing out some 
relevant directions for future research. 



II. MODEL 



We consider the following general form of the 2D s- 
DNLS equation in a rectangular geometry [TO] : 



. dU n ,n 
dz 
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\U n , r . 
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where z corresponds to the normalized propagation co- 
ordinate, 7 to an effective nonlinear parameter, and 
U n , m represents the (complex) electric field amplitude 
at site {n, m} in an array of N x M sites. A" m U nim = 

(frn+l,m + frn-l,m) + «(^n,m+l + ^n,m-l) defines the 

linear interaction between nearest-neighbor waveguides, 
where a is an anisotropy parameter allowing us to study 
different 2D scenarios from a completely isotropic lat- 
tice (a = 1) to a completely decoupled set of ID arrays 
(a = 0). Model possesses two conserved quantities, 
the power (norm) defined as 



i 2 

m | i 



(2) 



and the Hamiltonian (energy) defined as 

H = ~ ^2 [( U n+l,m + a#n,m+l) U*^ 
nm 

-|ln(l + |C/„, m | 2 ) + c.c. 
A. Stationary solutions 



(3) 



Stationary solutions to Eq. ([!]) are of the form 
U n: m(z) — u ni m exp (i\z), where u nj7n is a ^-independent 
(generally complex) amplitude and A corresponds to the 
propagation constant or frequency [10]. Extended sta- 
tionary solutions with constant amplitude (plane waves) 
exist in frequency bands, whose edges depend on the am- 
plitude \u ny7n \. The linear band region is obtained for 
low-amplitude plane waves (|u n>m | — » 0), and it is easy 
to show that they will exist in the region A G [—7 — 2(1 + 
a), —7 + 2(1 + a)]. Moreover, for high-amplitude plane 
waves (\u ny7n \ — >- 00), the nonlinear term saturates and 
the corresponding interval is [—2(1 + a), 2(1 + a)]. 

The fundamental localized stationary solutions to Eq. 
([!]), defined as solutions with real amplitudes u n ^ m having 
a single maximum distributed on the sites in one unit-cell 
of the lattice, were described for the isotropic case in [TP] . 
In general, we may identify one-site (OO), two-site hor- 
izontal (EO), two-site vertical (OE), and four-site (EE) 
solutions. Figs. [TJa)-(d) show typical profiles for these 
types of excitations. As discussed in [10], these in-phase 
localized solutions bifurcate from the fundamental linear 
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Figure 1: (Color online) Examples of spatial profiles for (a) 
one-site (00), (b) two-site horizontal (EO), (c) two-site verti- 
cal (OE), (d) four-site solutions (EE), (e) IS1, (f) IS2, (g) IS3- 
vertical, and (h) diagonal solutions, respectively, illustrated 
for an isotropic lattice (a = 1). 



mode (upper band edge) in the small- amplitude limit, 
and merge into the corresponding high-amplitude plane- 
wave mode in the infinite-power limit. As a consequence, 
they exist in the region A G [—7 + 2(1 + a), 2(1 + a)], so 
that the size of the existence region will be just 7. [A 
two-site diagonal solution [22 [see Fig. [ljh)] could also 
be excited but only in frequency regimes where solutions 
are very localized with A < 0, requiring large values of 7 
(for a = 1, 7 > 8.2)]. In addition, in certain parameter 
regimes there are also asymmetric, intermediate station- 
ary solutions, associated with exchange of stability be- 
tween the symmetric ones [10]. We may identify three 
types of such solutions, termed henceforth intermediate 
1 (IS1), intermediate 2 (IS2), and intermediate 3 (IS3), 
respectively. IS1 and IS3 both connect a stable two-site 
mode, OE/EO, with another simultaneously stable solu- 
tion: the one-site mode 00 (IS1), or the four-site mode 
EE (IS3) [see Figs.[TJe) and (g)]. Thus, in these cases the 
unstable intermediate solutions act as carriers of instabil- 
ity between the corresponding fundamental modes. The 
unstable IS2 solution exists when the two-site solutions 
are stable. It connects the unstable one-site solution 00 
with the likewise unstable four-site solution EE, and sta- 
bilizes the latter mode when reaching it [see Fig. [TJf)]. 

Figs. [2] shows fundamental properties of different sta- 
tionary solutions for an isotropic lattice with 7 = 4 
(similar curves were plotted and discussed in [10 , but 
we here include some additional examples to fa cilitat e 
the comparison with the full PN surfaces in Sec. Ill A). 
Power versus frequency curve is shown in Fig. [2] (a) , 
where different solutions repeatedly cross each other in 
the whole range of parameters. In Fig. [2] (a)-inset we 
plot AP/P = (Pi — Poo) /Poo as a function of A, for 
i = 00,EO/OE,EE,ISl,IS2JS3. This figure shows 
clearer how solutions cross each other including many 
intermediate solutions appearing in the region A ~ 1.3. 
Moreover, we also observe these crossings when plotting 
AH = Hi - Hoo versus power [Fig. [| (b)]. (We plot 
AH instead of H because Hamiltonian differences be- 
tween stationary solutions of the same power are gener- 



40 



30 



20 



10 




(a) 



0.5 



1.5 




Figure 2: (Color online) Properties of fundamental stationary 
solutions for 7 = 4 and a = 1: (a) Power versus frequency; 
(a)-inset: AP/P versus A; (b) AH versus power; (c) Stability 
versus power. Different solutions are plotted with different 
line-types: one-site (thick solid blue), two-site (thick dotted 
blue), four-site (thick dashed blue), IS1 (thin solid black), IS2 
(thin dotted black), and IS3 (thin dashed black), respectively. 



ally quite small, which is indeed a favorable scenario for 
moving solutions) . Inset in Fig. [2] (b) shows a zoom in 
the region P ~ 9.8 where the Hamiltonian of different 
fundamental solutions matches at some points. This was 
initially interpreted as a vanishing PN barrier [5 , which 
is actually not the case due to the intermediate solu- 
tions appearing in-between, raising the effective energy 
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barrier [10] [see Fig. [2] (b)-inset when thick blue lines co- 
incide]. However, good mobility would be expected close 
to these regions but, still, it will be strongly determined 
by the specific kick or perturbation given to the solution 
in order to put it in movement. 

In order to study the stability of stationary local- 
ized solutions we implement a standard method [6j [10] , 
obtaining all the stable/unstable linear perturbation- 
modes and their corresponding instability measure g = 
maxj[^{co>j}], where ujj are the oscillation frequencies of 
the linear eigenmodes. Thus, if g = the solution is (at 
least marginally) stable and if g > it is unstable. Sta- 
bility versus power for the same solutions as in Figs. [2] 
(a), (b) is shown in Fig. [2] (c) . (Note that in the corre- 
sponding figures 2-3 of [10] . the quantity G = — g 2 was 
used as instability measure.) From this figure, we can see 
how unstable intermediate solutions IS1, IS2, IS3 appear 
when two or three solutions (regarding OE and EO as 
different solutions) are simultaneously stable. Note that 
we never observed regions for simultaneously stable one 
and four-site solutions for isotropic coupling. However, 



as will be discussed further in Sec. [TVj in the anisotropic 
case such regions exist. 



B. Constraint method 

The constraint method allows us to construct energy 
surfaces connecting stationary solutions for a given value 
of power. In that sense, it helps us to effectively predict 
and interpret the dynamics across the lattice. Critical 
points will represent stationary solutions, and a coherent 
movement across the lattice should transform one solu- 
tion to the other by keeping the power constant [T2] . 
The method was originally introduced by Aubry and 
Cretegny [13] and then implemented by Savin et al. [23] 
to study the mobility of kinks in nonlinear Klein- Gordon 
lattices (see also Ref. [24] for a related approach to 
analyze travelling breathers in ID oscillator chains in 
terms of an effective Hamiltonian.) Lately, it was nu- 
merically implemented to analyze surface states in one- 
dimensional semi-infinite systems [21 . By adiabatically 
changing the amplitude in one particular site, specifi- 
cally chosen as the one after the main peak (the peak 
is at n c and the constrained amplitude at n c + 1), the 
one- and the two-site solutions could be connected. End- 
ing the sweep when u Uc = w nc +i an d the center of mass 
of the constrained solution, X, has varied from n c to 
n c + 0.5, a one-dimensional energy surface, H vs X, can 
be sketched. Technically speaking, the method used in 
Ref. [21] consists on eliminating one equation from the 
Newton-Raphson problem, the one of the constrained 
amplitude which is not anymore an unknown variable. 
However, as the power is kept constant, an equation for 
P is added and the frequency A becomes a variable com- 
pleting the variable-equations set. 

In the present work, for the 2D lattice we implement a 
more sophisticated method where we explicitly vary the 



center of mass instead of the amplitude. Since we are in- 
terested in the energy landscape around the fundamental 
stationary solutions, we will assume also the amplitude 
u n ^ m of the constrained solutions to be real and positive 
on the constraint sites. Then, from the definition of the 
center of mass coordinates in the horizontal and vertical 
directions, 



X 



Enm n \ u n,m\ 2 



and Y 



(4) 

we can solve for any pair of amplitudes A = u rhA ^ rnA and 
B u„. c 



at places {n Al m A } and {n^,m^}: 
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with j = n, m. Thus, as before the actual constraints will 
be in the amplitudes A and 5, but now we can tune the 
center of mass as wished, from a given stationary solution 
towards any other. In general, the constrained solutions 
obtained in this way will of course not be stationary solu- 
tions of the full system. To identify any stationary solu- 
tion - including the intermediate ones - we check whether 
the value of A obtained from the constrained NR scheme 
coincides with the frequency of a hypothetical true sta- 
tionary solution to the full Eq. ([!]) with the computed 
amplitude profile. Furthermore, constrained solutions al- 
low us to calculate their Hamiltonian, and therefore to 
construct an effective energy landscape. 
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Figure 3: (Color online) Array scheme showing the constraint 
locations used to obtain the energy surfaces below. 

The choice for positions {n Al m A } and {n Bl m B } is ev- 
idently not unique, but in order to most efficiently trace 
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out a smooth energy landscape connecting the funda- 
mental solutions (if it exists), the constraint sites should 
preferably be chosen within the unit cell where the am- 
plitudes are large. It turned out that, starting from a 
stationary one-site solution centered at {n c , m c }, in most 
cases the best option is to choose the first constraint at 
site {n c + l,m c } or {n c ,m c + 1}, and the other one at 
site {n c + 1, m c + 1}, as sketched in Fig. [3j The results 
shown in the following sections are obtained using these 
constraint sites. We also tried using constraint sites at 
{n c +l, m c } and {n c , m c +l} (i.e., along a diagonal); how- 
ever, with this choice we typically were not able to find 
the four-site EE solution when starting from the one-site 
00. Instead, the NR method with constraints on these 
sites generally converges, at X = n c + |, Y = m c + |, to 
a two-site diagonal solution (which as remarked in Sec. 
|II A| is an exact solution to the unconstrained equation 
only for large 7), which is not expected to be relevant for 
the mobility process. 



III. ISOTROPIC LATTICE 

A. Description of PN surfaces 

By exploring the parameter-space (7, P), we have 
identified regions where we were able to compute com- 
plete two-dimensional energy surfaces (using constraint 
as shown in Fig.[3| and regions where we cannot. In gen- 
eral, complete surfaces could not be obtained for large 
values of the nonlinearity constant, 7 > 6, when solu- 
tions get strongly localized. In some cases we could trace 
some specific mobility directions (e.g., along lattice axes 
or diagonal) connecting stationary solutions with almost 
equal Hamiltonian, but not the full two-dimensional sce- 
nario involving all fundamental solutions. Note also that, 
as remarked in [10], for large 7 and small power the s- 
DNLS model essentially behaves as a cubic DNLS model 
with effective nonlinearity 7P. Consequently, we could 
not find complete energy surfaces with all fundamental 
solutions for the cubic DNLS either (although surfaces 
involving the two-site diagonal solution in place of the 
EE solution could be obtained as mentioned in Sec. |IIB| 
for the cubic DNLS this two-site solution does exist as 
a true stationary solution above some threshold level of 
power [22]). 

On the other hand, for very small values of 7, the com- 
putation of energy surfaces becomes difficult for techni- 
cal reasons: due to the widening of the solutions, con- 
siderably larger lattices are needed to remove the influ- 
ence from boundary effects. Therefore, in the follow- 
ing, we will present the main phenomenology for the 
isotropic case (a = 1) found for intermediate values of 
7 7 3 < 7 < 5, where complete surfaces involving all 
four fundamental solutions are obtained for all values of 
the power. We will show results for the particular value 
7 = 4, but the scenario is found to be qualitatively the 
same for all 7 in this interval. 



From Fig. [2](c), we may identify essentially five differ- 
ent regimes where energy surfaces of qualitatively differ- 
ent nature should be expected, depending on the level of 
power. The first one is for low power, where (similarly 
to the cubic DNLS) the one-site solution is always sta- 
ble and the other fundamental stationary solutions are all 
unstable. The corresponding energy surface is illustrated 
in Fig. [4] (a), where the one-site solution yields the en- 
ergy minimum, the two-site solutions saddle points and 
the four-site solution the maximum. Note that, in this 
low-power regime, the surfaces for this value of 7 are still 
rather flat, and therefore some mobility may result if the 
one-site solution is kicked to overcome the barriers, in the 
axial as well as in the diagonal directions, as illustrated 
by Figs. 4 (a) and (c) in Ref. [10 j. (Anoth er exa mple of 
mobility in this regime is discussed in Sec. |IIIB below.) 



The saturable nature of the system becomes evident at 
higher powers. In the second regime, appearing for the 
first time when 9.27 < P < 9.55, the one- and the two- 
site solutions are stable simultaneously and, as shown in 
Fig. [4] (b) , these three points all correspond to local min- 
ima of the surface. Intermediate solutions IS1 connecting 
the one- with the two-site solutions in the horizontal and 
vertical directions appear as saddle points [white dots in 
Fig. [4] (b)]. Note that the energy landscape for the pa- 
rameter values in Fig. [4] (b) is almost flat between the 
one-site and two-site solutions in the axial directions, 
leading to the very good axial mobility shown in Fig. 
4 (d) of Ref. [10], while the maximum corresponding to 
the unstable four-site solution creates a too large effective 
barrier to overcome in the diagonal direction. 

The third power region is the one in which only the 
two-site solutions are stable. It appears for the first time 
when 9.55 < P < 10.04, and is illustrated in Fig. [| (c). 
Here, the stable two-site solutions correspond to two lo- 
cal minima of the surface, and the unstable one- and 
four-site solutions both to local maxima. The two unsta- 
ble solutions are connected by the intermediate solution 
IS2 (white dot at the surface), corresponding to a saddle 
point which for symmetry reasons (for the isotropic lat- 
tice) will lie along the diagonal connecting the unstable 
solutions. As will be illustrated in Sec. IIIB[ the easiest 
mobility in this case is expected to occur in a diagonal 
direction, connecting the two stable stationary solutions. 

The fourth regime corresponds to simultaneously sta- 
ble two- and four-site solutions. Such a regime appears 
for the first time when 10.04 < P < 10.32. If we 
further increase the power, we will find a similar re- 
gion with these three simultaneously stable solutions for 
34.5 < P < 36.25. In this regime, as illustrated in Fig. 
[4] (e), the two- and four-site solutions all correspond to 
minima, the unstable intermediate solutions IS3 to saddle 
points, and the one-site solution to a maximum. Thus, 
comparing Figs. [4] (b) and (e), we see that the structure 
of the potential has been completely inverted. Now, a 
big "hill" is located at the one-site solution, and as a 
consequence no mobility is expected involving this solu- 
tion. As will be shown below in Sec. IIIB[ the simplest 
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Figure 4: (Color online) Energy surfaces for 7 = 4, a = 1, in the five different power regimes discussed in the text : (a) 
P = 5 ; (b) P = 9.45; (c) P = 10; (d) P = 12; (e) P = 35.5. The center of mass {X, Y} for the four stationary solutions are: 
{8,8} (00), {8.5,8} (EO), {8,8.5} (OE), and {8.5,8.5} (EE). White dots denote local extrema corresponding to intermediate 
solutions. (System size N = M = 15 and fixed boundary conditions were used.) 
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mobility scenario will now be the one in which the two- 
site solutions travel through the four-site one across the 
lattice. 

The fifth regime, with the four-site solution being the 
only stable one, occurs for the first time when 10.32 < 
P < 34.5. As illustrated in Fig. |] (d), the four-site 
solution is now a minimum of the PN potential, while 
the other three unstable solutions correspond to saddle 
points (EO and OE) and maximum (00), respectively. 
(There are no intermediate solutions in this regime.) 
Thus, by increasing the power we have now reverted the 
surface compared to the low-power regime in Fig. [4] (a). 

Further increasing of power shows for 36.25 < P < 
43.6 a new region of stable two-site solutions, corre- 
sponding to the third regime. For 43.6 < P < 44.2 
the one- and two-site solutions are simultaneously sta- 
ble, so the scenario is equivalent to the second regime, 
followed again for P > 44.2 by the first regime, respec- 
tively. This second complete inversion thus corresponds 
to a regain of the low-power characteristics of the surface. 
Furthermore, repetition of these scenarios can be found 
for P > 118, and we confirmed, e.g., the existence of a 
complete, smooth surface, analogous to the one shown in 
Fig. [4] (d) for the fifth regime, for P = 120. 

B. Mobility dynamics for isotropic lattice 

To explicitly show the connection between the differ- 
ent types of energy surfaces described in Sec. |III A| and 
the mobility of localized solutions, we numerically inte- 
grate model by taking as initial condition a station- 
ary solution perturbed with a small kick: U niTn (0) = 
u ni7n ex.-p[ik x (n — n c ) + ik y (m — ra c )], where k x and k y 
correspond to the kick strength in the horizontal and 
vertical directions, respectively. If the surfaces would be 
completely flat, solutions should move even with an in- 
finitesimally small kick. However, as we have seen, in 
general surfaces are not flat due to the discreteness and 
the self-induced PN potential, and although the PN bar- 
rier can be very small in certain directions, it is generally 
non-zero. Therefore, in order to put a localized solution 
in movement, some amount of kinetic energy (represented 
by this kick) should be given to effectively overcome the 
energy barriers. 

We first discuss an example from the first, low-power, 
regime corresponding to the energy surface in Fig. [4] (a) 
(as mentioned above, different mobility examples from 
this regime were shown in Figs. 4 (a) and (c) in Ref. [10 ). 
Here, we took the (unstable) EO solution and kicked 
it with a very small value in both directions, yielding 
the dynamics numerically observed in Fig. [5] (a). In the 
^-direction, the kinetic energy is not sufficient to over- 
come the barrier created by the four-site solution, while 
in the x-direction it can move towards the minimum cor- 
responding to the one-site solution. The resulting dy- 
namics for the center-of-mass positions shows how the 
movement in the horizontal direction gets combined with 



oscillations in the vertical one. The dependence X vs z 
also clearly shows how the solution feels the potential in 
terms of its velocity: maximum velocities occur around 
integer lattices sites (where the potential is a minimum) 
while the minimum velocities occur close to middle points 
[see full line in Fig. [5] (a); compare also with the analo- 
gous scenarios for the kicked OO mode in Figs. 4 (a) and 
(c) of Ref. [10]]. 

In the second power regime, with energy surfaces as in 
Fig. [4] (b) , good mobility can be expected only in axial di- 
rections, with effective energy barrier determined by the 
intermediate solution between the one- and the two-site 
solutions as discussed in Ref. [10]. The dynamics in this 
power regime, with a very small kick applied only in the 
axial direction to one of the stable stationary solutions, 
will be as already illustrated in Fig. 4 (d) of Ref. [TO] : 
the solution moves very slowly and adiabatically traces 
the shape of the potential with a minimal velocity at the 
places of the intermediate solutions. 

In the third power regime, with surfaces as in Fig. [4] 
(c), a very interesting kind of mobility is observed: a di- 
agonal mobility between the (stable) horizontal and ver- 
tical two-site solutions as illustrated in Fig. [5] (b). The 
initial EO solution, kicked equally in the (—x)- and y- 
directions, gets sufficient kinetic energy to pass over the 
small barrier created by the intermediate IS2 solution. It 
then continues through the OE solution, passes another 
IS2 barrier, and then to the other EO solution shifted by 
one site in both directions. Although the potential con- 
necting these two solutions is not completely flat, there 
is a very good transport of energy in this direction, al- 
lowing mobility for more than one lattice diagonal in the 
considered 15 x 15 lattice. 

If we take a look at the surface in the fifth power regime 
[Fig. [4] (d)], we realize that an initial (unstable) one-site 
solution may move in any direction by slightly kicking it 
since it corresponds to a maximum. Fig. [5] (c) shows an 
example for kicking the one-site solution symmetrically in 
both directions in order to make it move passing through 
the four-site, i.e., a diagonal movement. The velocity has 
maximum in the minima of the potential (corresponding 
to the EE solution) and minimum in the potential max- 
ima (OO solution) (no intermediate solutions appear in 
this regime). 

Finally, in the fourth power regime, we see from the 
surface [Fig. [4] (e)] that the two- and four-site solutions 
are stable simultaneously, presenting an intermediate so- 
lution in-between that will define the effective energy bar- 
rier. This barrier is very small and, therefore, a very 
small kick is also needed. Figure [5] (d) shows the evo- 
lution starting from a two-site vertical solution, passing 
through the intermediate one and arriving to the four-site 
solution. 

As illustrated by the example in Fig. [5|a), the energy 
surfaces also provide nice intuitive interpretations to the 
dynamics of moving discrete solitons with additional per- 
turbations transverse to the direction of motion. In this 
example, the dynamics in the two orthogonal directions 
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Figure 5: (Color online) Examples of mobility dynamics in the propagation direction 2, corresponding to surfaces shown in Fig. 
[4] (a), (c), (d), and (e), respectively. In each subngure, top figures show profiles movement [(a) shows a colormap where colors 
were normalized to the maximum amplitude of each plot; this colormap also applies to (b)-(d)], and bottom figures center of 
mass evolution for X (full line) and Y (dashed line), (a) P = 5, k x = k y = 0.038; (b) P = 10, -k x = ky = 0.018; (c) P = 12, 
k x — k y — 0.015; (d) P — 35.5, k x — 0.015 and k y — 0. (Other parameter values are the same as in Fig. [4j) 



appear essentially independent of each other (propaga- 
tion in X while oscillating in Y). However, as we will 
illustrate with another example below, there are situa- 
tions where the particular surface topologies may lead 
to a more intricate interplay between the oscillatory and 
translational dynamics. 

We present in F ig. [6] a case with P = 44. With the 
notation from Sec. |III A[ this value of the power belongs 
to the second power regime, and the structure of the sur- 
face is phenomenologically identical to the one sketched 
in Fig. |4jb), but for a higher level of power. A picture 
of this energy landscape, periodically extended along the 
X-direction, is shown in Fig. [6] (a). In this figure, the 
darkest regions correspond to the positions for the min- 
ima (the 00 and EO/OE modes), while the brightest 
regions correspond to maxima (EE modes). Note also 
the positions of the saddle points (IS1, marked with ar- 
rows) in-between the 00 and EO/OE modes. 

The kick we choose for the vertical direction is very 
small, k y = 0.001, in order to disturb the horizontal 
movement only weakly. Since the potential shape in the 
F-direction - in a first approximation - will correspond 
to a harmonic potential, the profile will essentially per- 
form harmonic oscillations in the Indirection, with a fre- 
quency (determined by the curvature of the surface) that 
will vary only slightly as the mode translates in the X- 
direction. To look for possible resonances between the 
oscillatory and translational dynamics, we implement the 
following scheme: For fixed k x and k y we numerically in- 
tegrate model 0, starting with the 00 stationary solu- 
tion centered at position {8, 8}, and measure the z- values 
for which X is equal to 9, 10, 11 and 12, respectively. For 
these z- values we compute Y rn = Y(X = m), and plot in 
Fig.^b) as a function of the horizontal kick strength k x . 

The particular values of k x where symbols in Fig. [6^b) 
coincide for Y = 8 correspond to kick strengths where 
the natural frequency of the surface potential is com- 
mensurate with the translational velocity. In Fig. [6^c) 
we show some examples for these particular points. For 
k x = 0.0111 the solution makes two cycles before arriving 
to the next integer position in X. For k x = 0.0141 the 
solution makes one and a half cycles, for k x — 0.021 one 
cycle, for k x — 0.049 one half of cycle and for k x — 0.083 
one quarter of cycle, before arriving to the next integer 
position in X. As can be seen, around these points the 
oscillations apparently remain bounded and the dynam- 
ics stable, at least for long times. (Note that we do not 
see any visible phase-locking effects between the transla- 
tional and oscillatory motion here: the symbols in Fig. 
|6|b) seem to coincide only at isolated points and not in 
intervals. This may be an indication that the dissipative 
effects on the effective dynamics are very weak.) 

However, in other regimes we observe that the ampli- 
tude of the oscillations in the ^-direction rapidly grows. 
Two examples are illustrated in Fig.[6jd), corresponding 
to values of k x where the F-displacement in Fig.[6^b) has 
maxima. (Here, we have used a larger system of 41 x 41 
sites, in order to clearly identify the nature of this partic- 
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Figure 6: (Color online) Dynamics of an OO mode with 
P = 44 (7 = 4, a = 1), kicked with k y = 0.001 and vary- 
ing k x . (a) Extended energy surface. Darker (brighter) 
regions correspond to lower (higher) values of the energy 
H. (b) I9, Y10, Y11, Y12 versus k x represented by (black) 
filled circles, (grey) squares, (blue) diamonds and (red) 
triangles, respectively. (c) Y(z) versus X(z) for k x = 
0.0111, 0.0141, 0.021, 0.049 and 0.083 represented by thick 
(blue) solid, thin (black) dashed, (blue) dotted, thin (black) 
solid and thick (blue) dashed lines, respectively, (d) Dashed 
(blue) and solid (black) lines correspond to k x = 0.018 and 
k x = 0.035, respectively, for N = M = 41. 
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ular dynamics.) The initial increase of the ^-amplitude 
suggests a kind of unstable dynamics. However, as is 
clear from the longer simulation in Fig. [6jd), the oscilla- 
tion amplitude in the Y-direction remains bounded and 
thus the trajectory does not escape in the vertical direc- 
tion but remains trapped by the barrier. The typical dy- 
namics, thus, exhibits amplitude-modulated oscillations 
(beatings), with the largest amplitudes corresponding to 
maxima close to integer X where the surfaces are flat- 
ter. Thus, this kind of mobility is a signature of the 2D 
topology of the potential and it also validates the sur- 
faces computed with our method. (We performed several 
computations in different {P, /^j-regimes finding similar 
results as the ones presented in Fig. [6| . 

IV. EFFECTS OF WEAK ANISOTROPY 
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Figure 7: (Color online) Density plots of g versus P and a for 
7 = 4, for the (a) OO; (b) EO; (c) OE; and (d) EE solutions, 
respectively. Black means g — and lighter colors imply an 
increasing instability. 

Considering weak anisotropy by putting a < 1 in Eq. 
([!]) (implying a stronger coupling in the horizontal than 
in the vertical direction), the main qualitative modifica- 
tions as concerns the stationary solutions are shifts of 
the stability regions for the different solution types (ev- 
idently, the two-site horizontal, EO, and vertical, OE, 
solutions are now non-equivalent). Focusing our discus- 
sion on the parameter value 7 = 4 analyzed in detail 
for the isotropic case in the previous sections, the re- 
sults for the stability analysis are presented in Fig. [7[ 




Figure 8: (Color online) Energy surfaces for anisotropic lat- 
tices with 7 = 4: (a) a = 0.72, P = 19.5; (b) a = 0.7, 
P 20. 



We find that, for fixed a < 1 and increasing power, the 
first stability inversion appears between the one-site OO 
and the horizontal two-site EO solutions, and it occurs 
for lower values of the power than in the isotropic case. 
The first interval of stability of the four-site EE solution 
gets narrower when a decreases, mainly because its lower 
limit increases (with a corresponding increase of the up- 
per stability limit for the EO solution). The first stability 
region of the vertical two-site OE solution disappears for 
a < 0.96. On the other hand, its second stability regime 
for higher powers is enhanced, whereas the correspond- 
ing stability regime for its horizontal counterpart gets 
narrower and disappears for a < 0.89. 

With stronger anisotropy, some qualitatively new 
regimes not present in the isotropic case will appear. For 
a < 0.76, a new region of stability of the OO mode ap- 
pears around P ~ 20, leading to a small window of mul- 
tist ability where the OO, EO and EE modes are stable 
simultaneously. An example of an energy surface from 
this regime is shown in Fig. [8] (a), where also two inter- 
mediate solutions [white dots] can be found on the edges. 
There is also a regime with bistability of the OO and EE 
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modes only, as illustrated by the pseudo-potential surface 
in Fig. [8] (b). Note that, in this case, the intermediate 
solution does not leave the edge, in contrast to the IS2 
mode existing on the diagonal between the one- and four- 
site solutions for the isotropic case of simultaneously sta- 
ble EO/OE solution [Fig. [3] (c)]. Since anisotropy de- 
stroys the symmetry between the OE and EO solutions, 
bistability of both 2-site solutions disappears already for 
a small amount of anisotropy, and for smaller a the IS 
position can only be found on an edge of the potential 
surface. 

Thus, since no additional "valleys" were found to de- 
velop in the PN surfaces (all minima appear on the 
edges), the best mobility for anisotropic lattices should 
generally be expected to occur along lattice directions 
(at least in cases where diagonal coupling terms can be 
neglected, as assumed throughout this work). 

V. CONCLUSIONS 

In conclusion, we have studied in a deeper way 
the problem of mobility of localized modes in two- 
dimensional saturable discrete systems. We numerically 
implemented a constrained Newton-Raphson method to 
construct full Peierls-Nabarro energy surfaces, which ap- 
peared as very useful tools for predicting the dynamical 
properties of localized excitations. Although these sur- 
faces were never found to be completely flat (and there- 
fore the corresponding Peierls-Nabarro barriers is strictly 
never zero), parameter regimes and directions of good 
mobility were seen to correspond to smooth, flat parts of 
the surfaces. 

For the isotropic saturable model, five different surface 
topologies could be identified in different power regimes, 
depending on the stability properties of the different 
fundamental stationary solitons. By numerically study- 
ing the dynamics of perturbed stationary solutions, we 
showed how these different topologies yielded qualita- 
tively different kinds of optimal mobility, generally in 
axial or diagonal directions and with velocities varying 
according to the shape of the potential while the profile 
is propagating across the lattice. The energy surfaces 
were also found to be very useful for interpreting the dy- 
namics resulting from the interplay between translational 
and oscillatory motion in orthogonal directions. 



We also studied the effect of weak anisotropy, where 
the breaking of the symmetry of the system also yields 
energy surfaces with lower symmetry. In addition, we 
showed how this symmetry breaking could lead to a 
change of stability properties for the fundamental sta- 
tionary solutions, and as a result two new types of sur- 
face topologies were identified. From the shape of these 
surfaces, we predicted that the best mobility in the 
anisotropic case should generally occur along lattice di- 
rections. 

Thus, the method developed here appears to be a very 
good tool for the understanding of the mobility dynam- 
ics of localized excitations in higher-dimensional lattices, 
clearly showing the "real" energy barriers that localized 
solutions experience. We hope that this approach will be 
useful also for other types of nonlinar lattice models. 

Finally, from an application point of view, we have 
given further examples on how the saturable nature of 
nonlinearity promotes the generation of many different 
stability and mobility scenarios. The example presented 
in section [ill B| shows the possibility to steer the dynam- 
ics of the solution by using the natural frequencies of its 
own self-induced potential. As an application of this, we 
propose a discriminative optical switch which, for exam- 
ple, will only be activated when X and Y positions are 
- simultaneously - an integer number. This is one more 
example of nonlinear lattices being the key to success for 
controlling the propagation of light in future all-optical 
technologies. 
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